

## save stata data file (Rdata file not available at time) to working directory.  Here we preserved the name 576555173Canada LAPOP AmericasBarometer 2017 V1.0_W.dta

#read file
install.packages('readstata13')
library('readstata13')
abcan<-read.dta13("576555173Canada LAPOP AmericasBarometer 2017 V1.0_W.dta", convert.factors=F)

#demographics

abcan$sex<-abcan$q1
abcan$sex[abcan$sex==1]<-0
abcan$sex[abcan$sex==2]<-1
abcan$sex[abcan$sex==3]<-NA


abcan$ageyears<-abcan$q2


abcan$hhincome<-(abcan$q10new-1)/15



abcan$collegedegree[abcan$ed_can %in% c(7, 8)]<-1
abcan$collegedegree[!(abcan$ed_can %in% c(7, 8))]<-0



#openness to authoritarian governance items

abcan$limitopp<-(abcan$pop101-1)/6

abcan$churchillrev<-1-((abcan$ing4-1)/6)

abcan$democnotpref<-abcan$dem2
abcan$democnotpref[abcan$democnotpref %in% c(1, 3)]<-1
abcan$democnotpref[abcan$democnotpref==2]<-0

#half-sample only
abcan$milcoupcrime[abcan$jc10==1]<-1
abcan$milcoupcrime[abcan$jc10==2]<-0

#half-sample only
abcan$milcoupcorruption[abcan$jc13==1]<-1
abcan$milcoupcorruption[abcan$jc13==2]<-0

abcan$coup1<-abcan$milcoupcorruption
abcan$coup1[is.na(abcan$coup1)]<-0

abcan$coup2<-abcan$milcoupcrime
abcan$coup2[is.na(abcan$coup2)]<-0


abcan$coup<-abcan$coup1+abcan$coup2

abcan$closeparliament[abcan$jc15a==1]<-1
abcan$closeparliament[abcan$jc15a==2]<-0

abcan$disallowvote<-1-(abcan$d1-1)/9

abcan$disallowdemonstration<-1-(abcan$d2-1)/9

abcan$disallowrunforoffice<-1-(abcan$d3-1)/9

abcan$disallowspeeches<-1-(abcan$d4-1)/9

abcan$toomuchpressfreedombinary[abcan$lib1==3]<-1
abcan$toomuchpressfreedombinary[abcan$lib1 %in% c(1, 2)]<-0

abcan$toomuchexpressionfreedombinary[abcan$lib2b==3]<-1
abcan$toomuchexpressionfreedombinary[abcan$lib2b %in% c(1, 2)]<-0

abcan$toomuchpolitfreedombinary[abcan$lib2c==3]<-1
abcan$toomuchpolitfreedombinary[abcan$lib2c %in% c(1, 2)]<-0


#left economic attitude item


abcan$incomeineq<-(abcan$ros4-1)/6

#cultural conservatism

#nones weren't asked, so are coded 0 for attendance
abcan$religattend<-1-((abcan$q5a-1)/4)
abcan$religattend[is.na(abcan$religattend)]<-0

abcan$religimport<-1-((abcan$q5b-1)/3)

abcan$againstsamesexmarriage<-1-((abcan$d6-1)/9)

abcan$punishcriminals<-(abcan$aoj22new-1)/6


#control variables


abcan$trustparliament<-(abcan$b13-1)/6

abcan$trustparties<-(abcan$b21-1)/6

abcan$trustlocalgov<-(abcan$b32-1)/6

abcan$trustelections<-(abcan$b47a-1)/6



abcan$satisdemoc<-1-((abcan$pn4-1)/3)


abcan$polint<-1-((abcan$pol1-1)/3)

abcan$follownews<-1-((abcan$gi0-1)/4)

#descriptives

wtd.mean(abcan$coup, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$coup, weights=abcan$wt, na.rm=T))


wtd.mean(abcan$closeparliament, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$closeparliament, weights=abcan$wt, na.rm=T))

wtd.mean(abcan$churchillrev, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$churchillrev, weights=abcan$wt, na.rm=T))


abcan$intol<-with(abcan, rowMeans(data.frame(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches)), na.rm=T)

wtd.mean(abcan$intol, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$intol, weights=abcan$wt, na.rm=T))

abcanintolitems<-subset(abcan, select=c(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches))

psych::alpha(abcanintolitems)




abcan$toomuchfreedom<-with(abcan, rowMeans(data.frame(toomuchpressfreedombinary, toomuchexpressionfreedombinary, toomuchpolitfreedombinary)),na.rm=T)

wtd.mean(abcan$toomuchfreedom, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$toomuchfreedom, weights=abcan$wt, na.rm=T))

abcantoomuchfreedomitems<-subset(abcan, select=c(toomuchpressfreedombinary, toomuchexpressionfreedombinary, toomuchpolitfreedombinary))

psych::alpha(abcantoomuchfreedomitems)


wtd.mean(abcan$limitopp, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$limitopp, weights=abcan$wt, na.rm=T))


wtd.mean(abcan$democnotpref, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$democnotpref, weights=abcan$wt, na.rm=T))


abcan$authgov<-  with(abcan, rowMeans(data.frame(coup, closeparliament, churchillrev, intol, toomuchfreedom, limitopp, democnotpref)),na.rm=T)

wtd.mean(abcan$authgov, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$authgov, weights=abcan$wt, na.rm=T))


abcanauthgovitems<-subset(abcan, select=c(coup, closeparliament, churchillrev, intol, toomuchfreedom, limitopp, democnotpref))

psych::alpha(abcanauthgovitems)




abcan$cultright<-  with(abcan, rowMeans(data.frame(religattend, religimport, againstsamesexmarriage, punishcriminals)),na.rm=T)


abcan$cultrightnoattend<-  with(abcan, rowMeans(data.frame(religimport, againstsamesexmarriage, punishcriminals)),na.rm=T)

wtd.mean(abcan$cultright, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$cultright, weights=abcan$wt, na.rm=T))

abcancultrightitems<-subset(abcan, select=c(religattend, religimport, againstsamesexmarriage, punishcriminals))

psych::alpha(abcancultrightitems)



abcan$econleft<-abcan$incomeineq

wtd.mean(abcan$econleft, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$econleft, weights=abcan$wt, na.rm=T))


abcan$confinstit<-with(abcan, rowMeans(data.frame(trustelections, trustparties, trustparliament, trustlocalgov)),na.rm=T)

wtd.mean(abcan$confinstit, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$confinstit, weights=abcan$wt, na.rm=T))

abcanconfinstititems<-subset(abcan, select=c(trustelections, trustparties, trustparliament, trustlocalgov))

psych::alpha(abcanconfinstititems)


abcan$poleng<-with(abcan, rowMeans(data.frame(polint, follownews)),na.rm=T)

wtd.mean(abcan$poleng, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$poleng, weights=abcan$wt, na.rm=T))

abcanpolengitems<-subset(abcan, select=c(polint, follownews))

psych::alpha(abcanpolengitems)



wtd.mean(abcan$satisdemoc, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$satisdemoc, weights=abcan$wt, na.rm=T))



wtd.mean(abcan$collegedegree, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$collegedegree, weights=abcan$wt, na.rm=T))



wtd.mean(abcan$ageyears, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$ageyears, weights=abcan$wt, na.rm=T))

wtd.mean(abcan$sex, weights=abcan$wt, na.rm=T)

sqrt(wtd.var(abcan$sex, weights=abcan$wt, na.rm=T))




#main analyses


library(mice)
library(miceadds)
library(survey)

xtract.svymi <- function(z){
  #require(norm)
  m2 <- length(z)
  ns <- nobs(z[[1]]) #first dataset
  betas <- lapply(z, FUN = function(x){coef(x)} )
  sess <- lapply(z, FUN = function(x){vcov(x)}) #the vcov
  #ses <- lapply(sess, FUN = function(x){sqrt(diag(x))}) #the se, not the complete vcov
  
  dm <- lapply(z, FUN = function(x){model.matrix(x)})
  r.1 <- sapply(1:m2, function(row) cor(dm[[row]] %*% betas[[row]], z[[row]]$y))
  fisher.r2z <- function(r) {0.5 * (log(1+r) - log(1-r))}
  zz <- mean(fisher.r2z(r.1))
  r2 <- psych::fisherz2r(zz)^2 #in psych
  R2 <- mean(r2)
  
  test1 <- summary(pool_mi(betas, sess))
  
  test1$term <- row.names(test1)
  test2 <- test1[,c(8,1,2,4)] #4 is the pvalue
  
  row.names(test2) <- c()
  names(test2) <- c('term','estimate','std.error', 'p.value')
  
  tr <- createTexreg(
    coef.names = test2$term, 
    coef = test2$estimate, 
    se = test2$std.error, 
    pvalues = test2$p.value,
    gof.names = c("R2", "Nobs"), 
    gof = c(R2, ns),
    gof.decimal = c(T,F)
  )
  
}


abcan<-subset(abcan, select=c(sex, ageyears, collegedegree, hhincome, closeparliament, coup, churchillrev, disallowvote, 
                            disallowdemonstration, disallowrunforoffice, 
                            disallowspeeches, toomuchpressfreedombinary, toomuchexpressionfreedombinary, 
                            toomuchpolitfreedombinary, limitopp, democnotpref, incomeineq, religattend, religimport, againstsamesexmarriage, punishcriminals, 
                            trustelections, trustparties, trustparliament, trustlocalgov, polint, follownews, satisdemoc, estratopri, wt, upm))



abcani<-mice(abcan, m=20)

longabcan <- complete(abcani, action='long', include=TRUE)

longabcan$intol<-with(longabcan, rowMeans(data.frame(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches)))

longabcan$toomuchfreedom<-with(longabcan, rowMeans(data.frame(toomuchpressfreedombinary, toomuchexpressionfreedombinary, toomuchpolitfreedombinary)))

longabcan$authgov<-  with(longabcan, rowMeans(data.frame(coup, closeparliament, churchillrev, intol, toomuchfreedom, limitopp, democnotpref)))


longabcan$cultright<-  with(longabcan, rowMeans(data.frame(religattend, religimport, againstsamesexmarriage, punishcriminals)))


longabcan$econleft<-  longabcan$incomeineq

longabcan$age1s<-with(longabcan, (ageyears-min(ageyears, na.rm=T))/(max(ageyears, na.rm=T)-min(ageyears, na.rm=T)) )

longabcan$confinstit<-with(longabcan, rowMeans(data.frame(trustelections, trustparties, trustparliament, trustlocalgov)))

longabcan$poleng<-with(longabcan, rowMeans(data.frame(polint, follownews)))


longabcan$cultrightplussd<-with(longabcan, cultright-(mean(cultright, na.rm=T)+sd(cultright, na.rm=T)))

longabcan$cultrightminussd<-with(longabcan, cultright-(mean(cultright, na.rm=T)-sd(cultright, na.rm=T)))

longabcan$econleftplussd<-with(longabcan, econleft-(mean(econleft, na.rm=T)+sd(econleft, na.rm=T)))

longabcan$econleftminussd<-with(longabcan, econleft-(mean(econleft, na.rm=T)-sd(econleft, na.rm=T)))

longabcan$cultrightscl <- with(longabcan, (cultright-mean(cultright, na.rm=T))/(2*sd(cultright, na.rm=T)))

longabcan$econleftscl <- with(longabcan, (econleft-mean(econleft, na.rm=T))/(2*sd(econleft, na.rm=T)))

longabcan$polengscl <-  with(longabcan, (poleng-mean(poleng, na.rm=T))/(2*sd(poleng, na.rm=T)))


longabcan$age1sscl <- with(longabcan, (age1s-mean(age1s, na.rm=T))/(2*sd(age1s, na.rm=T)))


longabcan$satisdemocscl <- with(longabcan, (satisdemoc-mean(satisdemoc, na.rm=T))/(2*sd(satisdemoc, na.rm=T)))


longabcan$confinstitscl <- with(longabcan, (confinstit-mean(confinstit, na.rm=T))/(2*sd(confinstit, na.rm=T)))


longabcan$hhincomescl <- with(longabcan, (hhincome-mean(hhincome, na.rm=T))/(2*sd(hhincome, na.rm=T)))


longabcan$econextrem<-abs(with(longabcan, (econleft-mean(econleft, na.rm=T))))

longabcan$econextremscl<-with(longabcan, (econextrem-mean(econextrem, na.rm=T))/(2*sd(econextrem, na.rm=T)))

abcani2 <- as.mids(longabcan)



abcan_imp_list <- imputationList(lapply(1:abcani2$m, function(n) mice::complete(abcani2, action = n)))

abcandesign<-svydesign(~upm,strata=~estratopri, nest=TRUE, fpc=NULL, data=abcan_imp_list, weights=~wt)


abcanreg1<- with(abcandesign,svyglm(authgov~cultrightscl+econleftscl))


abcanreg2<- with(abcandesign,svyglm(authgov~cultrightscl+econleftscl+cultrightscl:econleftscl))




abcanreg3<- with(abcandesign,svyglm(authgov~cultrightscl+econleftscl+sex+age1sscl+collegedegree+hhincomescl+confinstitscl+polengscl+satisdemocscl))


abcanreg4<- with(abcandesign,svyglm(authgov~cultrightscl+econleftscl+sex+age1sscl+collegedegree+hhincomescl+confinstitscl+polengscl+satisdemocscl+cultrightscl:econleftscl))




abcanregright<- with(abcandesign,svyglm(authgov~cultrightplussd+econleftminussd+cultrightplussd:econleftminussd))


abcanregleft<- with(abcandesign,svyglm(authgov~cultrightminussd+econleftplussd+cultrightminussd:econleftplussd))


abcanregprot<- with(abcandesign,svyglm(authgov~cultrightplussd+econleftplussd+cultrightplussd:econleftplussd))


abcanregfree<- with(abcandesign,svyglm(authgov~cultrightminussd+econleftminussd+cultrightminussd:econleftminussd))

abcanout1 <- xtract.svymi(abcanreg1)
abcanout2 <- xtract.svymi(abcanreg2)
abcanout3 <- xtract.svymi(abcanreg3)
abcanout4 <- xtract.svymi(abcanreg4)


htmlreg(list(abcanout1,abcanout2,abcanout3,abcanout4), digits = 3,file = "lapopCAN17.html",custom.coef.names = c("Constant","Cultural Conservatism","Left Economic Attitudes","Cultural Conservatism x Left Economic Attitudes","Female","Age","College Degree","Household Income","Confidence in Institutions","Political Engagement","Satisfaction with Democracy"),caption.above = T,caption = "Table D-16: Predictors of Openness to Authoritarian Governance in LAPOP Canada-2017 Dataset",custom.note = ("%stars. <br> Continuous predictors are mean-centered and scaled by two standard deviations. Binary predictors are coded 0-1.<br>The outcome variable is coded to range from 0.0 to 1.0 <br> Missing data were replaced using multiple imputation with results pooled across 20 copies of the dataset."),custom.gof.names = c("R<sup>2</sup>","Observations"))

abcancultest <- summary(pool(abcanreg1),conf.int = T)[2,2]
abcancultlower <- summary(pool(abcanreg1),conf.int = T)[2,7]
abcancultupper <- summary(pool(abcanreg1),conf.int = T)[2,8]

abcaneconest <- summary(pool(abcanreg1),conf.int = T)[3,2]
abcaneconlower <- summary(pool(abcanreg1),conf.int = T)[3,7]
abcaneconupper <- summary(pool(abcanreg1),conf.int = T)[3,8]

abcanintest <- summary(pool(abcanreg2),conf.int = T)[4,2]
abcanintlower <- summary(pool(abcanreg2),conf.int = T)[4,7]
abcanintupper <- summary(pool(abcanreg2),conf.int = T)[4,8]


abcanrightest <- summary(pool(abcanregright),conf.int = T)[1,2]
abcanrightlower <- summary(pool(abcanregright),conf.int = T)[1,7]
abcanrightupper <- summary(pool(abcanregright),conf.int = T)[1,8]


abcanleftest <- summary(pool(abcanregleft),conf.int = T)[1,2]
abcanleftlower <- summary(pool(abcanregleft),conf.int = T)[1,7]
abcanleftupper <- summary(pool(abcanregleft),conf.int = T)[1,8]


abcanprotest <- summary(pool(abcanregprot),conf.int = T)[1,2]
abcanprotlower <- summary(pool(abcanregprot),conf.int = T)[1,7]
abcanprotupper <- summary(pool(abcanregprot),conf.int = T)[1,8]


abcanfreeest <- summary(pool(abcanregfree),conf.int = T)[1,2]
abcanfreelower <- summary(pool(abcanregfree),conf.int = T)[1,7]
abcanfreeupper <- summary(pool(abcanregfree),conf.int = T)[1,8]



abcancultb3<-extract.cis(abcanreg3)[2,1]
abcanpolengb3<-extract.cis(abcanreg3)[9,1]
abcanageb3<-extract.cis(abcanreg3)[5,1]
abcancollegeb3<-extract.cis(abcanreg3)[6,1]
abcanconfinstitb3<-extract.cis(abcanreg3)[8,1]


#additional analyses without imputation, scales with all avail items and listwise deletion in analysis.  Reported in Part E of SOM.


abcan$cultrightscl<-with(abcan, (cultright-mean(abcan$cultright, na.rm=T))/(2*sd(cultright, na.rm=T)))

abcan$econleftscl<-with(abcan, (econleft-mean(abcan$econleft, na.rm=T))/(2*sd(econleft, na.rm=T)))


abcan$cultrightplussd<-with(abcan, cultright-(mean(cultright, na.rm=T)+sd(cultright, na.rm=T)))

abcan$cultrightminussd<-with(abcan, cultright-(mean(cultright, na.rm=T)-sd(cultright, na.rm=T)))

abcan$econleftplussd<-with(abcan, econleft-(mean(econleft, na.rm=T)+sd(econleft, na.rm=T)))

abcan$econleftminussd<-with(abcan, econleft-(mean(econleft, na.rm=T)-sd(econleft, na.rm=T)))


abcandesign<-svydesign(~upm, strata=abcan$estratopri, nest=TRUE, fpc=NULL, data=abcan, weights=abcan$wt)


reg1abcanl<- svyglm(authgov~cultrightscl+econleftscl, design=abcandesign)


reg2abcanl<- svyglm(authgov~cultrightscl+econleftscl+cultrightscl:econleftscl, design=abcandesign)

abcanlregright<- svyglm(authgov~cultrightplussd+econleftminussd+cultrightplussd:econleftminussd, design=abcandesign)


abcanlregleft<- svyglm(authgov~cultrightminussd+econleftplussd+cultrightminussd:econleftplussd, design=abcandesign)


abcanlregprot<- svyglm(authgov~cultrightplussd+econleftplussd+cultrightplussd:econleftplussd, design=abcandesign)


abcanlregfree<- svyglm(authgov~cultrightminussd+econleftminussd+cultrightminussd:econleftminussd, design=abcandesign)



abcanlcultest<-summary(reg1abcanl)$coef[2,1]
abcanlcultlower<-confint(reg1abcanl)[2,1]
abcanlcultupper<-confint(reg1abcanl)[2,2]

abcanleconest<-summary(reg1abcanl)$coef[3,1]
abcanleconlower<-confint(reg1abcanl)[3,1]
abcanleconupper<-confint(reg1abcanl)[3,2]

abcanlintest<-summary(reg2abcanl)$coef[4,1]
abcanlintlower<-confint(reg2abcanl)[4,1]
abcanlintupper<-confint(reg2abcanl)[4,2]


abcanlrightest<-summary(abcanlregright)$coef[1,1]
abcanlrightlower<-confint(abcanlregright)[1,1]
abcanlrightupper<-confint(abcanlregright)[1,2]

abcanlleftest<-summary(abcanlregleft)$coef[1,1]
abcanlleftlower<-confint(abcanlregleft)[1,1]
abcanlleftupper<-confint(abcanlregleft)[1,2]

abcanlprotest<-summary(abcanlregprot)$coef[1,1]
abcanlprotlower<-confint(abcanlregprot)[1,1]
abcanlprotupper<-confint(abcanlregprot)[1,2]

abcanlfreeest<-summary(abcanlregfree)$coef[1,1]
abcanlfreelower<-confint(abcanlregfree)[1,1]
abcanlfreeupper<-confint(abcanlregfree)[1,2]

